R Notebook

What are the most frequent minor variants in the high-confidence samples?

source("./scripts/startup.R")

The following packages are a base install and will not be unloaded:



The following packages were not previously loaded:



Loading required package: pacman

patient_var_30 = read_feather("processing/patient_var_30.arrow")
maf_histogram = patient_var_30 %>% ggplot(aes(minorfreq)) + geom_histogram(binwidth = 0.01) + theme_pubr()
quantile(patient_var_30$minorfreq, probs = c(.25,.5,.75), na.rm = T)

       25%        50%        75% 
0.01252470 0.01717647 0.02788318 

ggsave("ggsave/maf_histogram.pdf", maf_histogram, width = 3, height = 3)

patient_var_30_for_rank = patient_var_30 %>% drop_na(gene) %>% filter(gene !="") %>%
  mutate(label = paste0(ref_sym, aapos), 
         refnt = str_sub(ref_codon, start = codon_pos, end = codon_pos))
bp = c("A", "T", "C", "G")
lineage_defining_mutations = fread("20220103-TRACE-LineageDefinitions-v9.1.txt", data.table = F) %>%
  filter(nt_ref %in% bp & 
           nt_alt %in% bp & 
           variation_type == "SNP") %>% select(nt_pos, nt_ref, nt_alt)
ldm_map = patient_var_30_for_rank %>% 
  mutate(ldm  = ifelse(ntpos %in% lineage_defining_mutations$nt_pos, TRUE, FALSE)) %>% 
  select(gene, label, ldm) %>% distinct %>% group_by(gene, label) %>% arrange(-ldm) %>%
  filter(row_number()==1)
  

highly_shared_sites_ranked = patient_var_30_for_rank %>% drop_na(gene) %>% filter(gene !="") %>%
  mutate(label = paste0(ref_sym, aapos), 
         refnt = str_sub(ref_codon, start = codon_pos, end = codon_pos)) %>%
  group_by(mcov_id, gene, label) %>%
  tally() %>% 
  arrange(-n) %>% ungroup %>% group_by(gene, label) %>% 
  summarize(label_count = n()) %>% ungroup() %>% arrange(-label_count) %>% mutate(rank = 1:nrow(.)) %>%
  left_join(ldm_map) %>% unique()

# highly_shared_sites_ranked_50 = highly_shared_sites_ranked %>% filter(rank <= 25)
# 
# ((
#   hss_maf_50<-patient_var_30 %>% filter(paste0(gene, aapos) %in% paste0(highly_shared_sites_ranked_50$gene, parse_number(highly_shared_sites_ranked_50$label))) %>% 
#     ggplot(aes(x = as.factor(ntpos), y = minorfreq)) + 
#     #geom_point(alpha = 0.2, position = position_jitter(width = 0.1)) + theme_bw() + 
#     geom_boxplot(outlier.shape = NA, alpha =0.5) + theme_bw() +
#     theme(axis.text.x = element_text(angle=90)) + #geom_boxplot(outlier.alpha = 0) + 
#     xlab("Nucleotide position") + ylab("Minor allele freq") ))
# 
# ggsave("ggsave/hss_maf_50.pdf", hss_maf_50, width = 5, height = 3)


#%>% 
 # left_join(patient_var_30 %>% select(ntpos, gene, ref_sym, aapos, ref_codon, codon_pos) %>% distinct)

# ((
#   highly_shared_sites<- patient_var_30 %>% group_by(ntpos) %>% tally() %>% arrange(desc(n)) %>% 
#        drop_na(ntpos) %>% slice_max(n, n = 35) %>% pull(ntpos) %>% as.numeric()
# ))


sample_n = length(unique(patient_var_30$mcov_id))
#for (this_gene in gene_limits$gene_id) {
  plot_ranks = highly_shared_sites_ranked %>% filter(rank <= 35) %>% 
    ggplot(aes(rank, label_count, label = label)) + 
  geom_line(data = highly_shared_sites_ranked %>% 
              select(!gene) %>% filter(rank <= 35), aes(rank, label_count), color = "grey") +
  scale_y_continuous(trans = "log2", breaks = 2^seq(0,12,1),
                     sec.axis = sec_axis(~./sample_n, labels = scales::label_percent(),
                                           breaks = 2^seq(0,8,1)/100)) + theme_pubr() +
  geom_point(aes(rank, label_count), color = "grey") +
  geom_text_repel(aes(label = label, color = ldm), ylim = c(6,11), xlim = c(NA, NA), angle = 90,
                  segment.size = 0.6, segment.curvature = -1e-20, segment.linetype = 3,
                  max.overlaps = Inf) + coord_cartesian(clip = "off") +
    scale_color_manual(values = c("salmon", "grey")) +
    facet_wrap(~gene)
  
  plot_ranks
#}
ggsave("ggsave/plot_ranks.pdf", plot_ranks, height = 5.5, width = 7)

# Across genome


lineage_defining_mutations

plot_data_LDM_prop = patient_var_30 %>% group_by(MCoVNumber) %>% drop_na(ntpos) %>% 
  summarize(iSNV_count = n(), iSNV = paste0(ntpos, collapse = ","), 
            iSNV_LDM_count = length(ntpos[ntpos %in% lineage_defining_mutations$nt_pos]),
            iSNV_LDM = paste0(ntpos[ntpos %in% lineage_defining_mutations$nt_pos],
                              collapse = ","), prop = iSNV_LDM_count/iSNV_count)

plot_data_LDM_prop %>% ggplot(aes(as.factor(iSNV_count), prop)) + geom_violin(draw_quantiles=c(0.5))

patient_var_30_reversions = patient_var_30 %>% mutate(reversion = minor == str_sub(ref_codon, 
                            start = codon_pos, end = codon_pos)) %>% drop_na(reversion)

# How many of our alleles were reversions?
patient_var_30_reversions %>% summarize(proportion_reversion = sum(reversion==T)/
                                          nrow(patient_var_30))

# How many of our alleles that are reversions are at LDM sites?
patient_var_30_reversions %>% mutate(ldm = ifelse(ntpos %in% lineage_defining_mutations$nt_pos,1,0)) %>% 
  summarize(proportion_reversion = sum(ldm == 1 & reversion == T)/sum(reversion==T))

# obscure the consensus in the VCF
problem_sites_global = fread("problematic_sites_sarsCov2_v8-20211027.vcf", skip = 88) %>%
  filter(FILTER != "mask") %>%
  filter(!grepl('single_src|nanopore', INFO)) %>% pull(POS)
problem_sites_houston = fread("problematic_sites_sarsCov2_v8-20211027.vcf", skip = 88) %>%
  filter(grepl('Houston', INFO)) %>% pull(POS)
problem_sites = unique(c(problem_sites_global, problem_sites_houston))
problem_sites

 [1]    76    78   320   538   660  1001  1814  1947  2087  5393  5498  6309  6310  6312  6866  7396
[17] 12685 13512 13513 13514 13686 13693 18505 18506 19338 19344 19369 19732 21302 21304 21305 21658
[33] 22393 22410 22488 22515 23144 24673 25798 26709 27792 28881 28882 28883 29378  8658 12698 15103
[49] 16130 16132 22892

# UNCOMMENT BELOW IF NECESSARY TO REGENERATE DATA
# consensus <-fread("consensus_minor_changes_20220713.csv", data.table = F) %>% filter(!ntpos %in% problem_sites) %>% filter(!ntpos %in% c(1:265)) %>% filter(!ntpos>29674)
# 
# 
# 
# consensus 
# 
# #
# consensus_count = consensus %>% 
#   mutate(MCoVNumber = mcov_reformat(name)) %>% 
#   select(MCoVNumber, major, ntpos, refnt) %>% 
#   filter(major != refnt) %>% 
#   unique() %>%
#   group_by(ntpos) %>%
#   summarize(MCoVNumber = unique(MCoVNumber)) 
# 
# runs = fread("sample_date_and_run.csv", data.table = F) %>%
#     mutate(MCoVNumber = mcov_reformat(mcov_id))
# 
# consensus_runs = consensus_count %>% left_join(runs) %>% drop_na() %>%
#   select(ntpos, MCoVNumber_possible_contaminant = MCoVNumber, run_consensus = run_group) #has each of the consensus mutations with the MCOV
# 
# # uncomment below if necessary
# #patient_var_30_contaminated = patient_var_30 %>% left_join(consensus_runs) # intensive
# #write_feather(patient_var_30_contaminated, "processing/patient_var_30_contaminated.arrow")
# 
patient_var_30_contaminated = read_feather("processing/patient_var_30_contaminated.arrow")
contaminated_df = patient_var_30_contaminated %>% mutate(ref_nt = str_sub(ref_codon,
                            start = codon_pos, end = codon_pos)) %>%
  select(MCoVNumber, run_group, ntpos, ref_nt, major, minor, MCoVNumber_possible_contaminant, run_consensus)

contaminated_df_tallied = contaminated_df %>%
  filter(run_group == run_consensus) %>%
  group_by(MCoVNumber) %>%
  summarize(n_var_contam = length(unique(ntpos)), MCoVNumber_possible_contaminant) %>%
  group_by(MCoVNumber, MCoVNumber_possible_contaminant) %>%
  summarize(n_var_contam, single_sample = n(), prop_single_sample = single_sample / n_var_contam)

`summarise()` has grouped output by 'MCoVNumber'. You can override using the `.groups` argument.`summarise()` has grouped output by 'MCoVNumber', 'MCoVNumber_possible_contaminant'. You can override using the `.groups` argument.

contaminated_df_tallied_top = contaminated_df_tallied %>% group_by(MCoVNumber) %>% top_n(n=1) %>%
  filter(row_number()==1)

Selecting by prop_single_sample

#Note prop_single_sample is single_sample / n_var_contam (proportion of number of minor variants that resememble consensus in another sample / number of minor variants that resemble any consensus in the run)
# prop single contam is the same as above but the denom is out of n_var
patient_counts_30 = read_feather("processing/patient_counts_30.arrow")
patient_counts_30_anno = patient_counts_30 %>% left_join(contaminated_df_tallied_top) %>%
  mutate(run_num = as.numeric(str_replace(run_group, "Run_", ""))) %>%
  replace(is.na(.),0)

Joining, by = "MCoVNumber"

# Plot: n_var x prob that it's from a single source
plot_contam = patient_counts_30_anno %>% mutate(prop_contam = n_var_contam/n_var, prop_single_contam = single_sample/n_var) 

plot_contam %>% select(n_var, admitted_hospital, prop_contam, prop_single_contam) %>% 
  pivot_longer(cols= c(prop_contam, prop_single_contam), names_to = "type", values_to = "prop")  %>%
  ggplot(aes(x = admitted_hospital, y = prop, 
             fill = type)) + 
  geom_violin(draw_quantiles = c(0.25,0.5,0.75))

median_contam_contam = quantile(plot_contam %>% filter(n_var > 0) %>% 
                                  pull(prop_contam), probs = c(0.5))
print(median_contam_contam)

50% 
0.5 

plot_contam_contam = plot_contam %>% 
  select(n_var, admitted_hospital, prop_contam, prop_single_contam) %>% 
  ggplot(aes(x = n_var, y = prop_contam, 
             fill = admitted_hospital, color = admitted_hospital)) + 
  geom_point(alpha = 0.02) + geom_smooth() + theme_pubr() +
  geom_hline(yintercept = median_contam_contam, linetype = "dashed")

plot_contam_contam = ggMarginal(plot_contam_contam, groupColour = T, groupFill = T, type = "violin",
           draw_quantiles = c(0.25, 0.5, 0.75))

`geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'Warning: Removed 450 rows containing non-finite values (`stat_smooth()`).`geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'Warning: Removed 450 rows containing non-finite values (`stat_smooth()`).Warning: Removed 450 rows containing missing values (`geom_point()`).

median_single_contam = quantile(plot_contam %>% filter(n_var > 0) %>% 
                                  pull(prop_single_contam), probs = c(0.5))
print(median_single_contam)

      50% 
0.3333333 

plot_contam_single_contam = plot_contam %>% 
  select(n_var, admitted_hospital, prop_contam, prop_single_contam) %>% 
  ggplot(aes(x = n_var, y = prop_single_contam, 
             fill = admitted_hospital, color = admitted_hospital)) + 
  geom_point(alpha = 0.02) + geom_smooth() + theme_pubr() + 
  geom_hline(yintercept = median_single_contam, linetype = "dashed")

plot_contam_single_contam = ggMarginal(plot_contam_single_contam, 
                                       groupColour = T, groupFill = T, 
                                       type = "violin",
           draw_quantiles = c(0.25, 0.5, 0.75))

`geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'Warning: Removed 450 rows containing non-finite values (`stat_smooth()`).`geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'Warning: Removed 450 rows containing non-finite values (`stat_smooth()`).Warning: Removed 450 rows containing missing values (`geom_point()`).

plot_contam_arranged = ggarrange(plot_contam_contam, plot_contam_single_contam,
                                 align = "v",
                                 labels = "AUTO")
plot_contam_arranged
ggsave("ggsave/plot_contam_arranged.pdf", plot_contam_arranged, height = 3, width = 6)

#consensus_sites_lowct %>% filter(refnt == major)
patient_counts_30 = read_feather("processing/patient_counts_30.arrow")


# criteria to count the mutation
# consensus = if ntpos matches in the list above, then 
# the minor will also match the nt_ref (i.e. reversion) or the nt_alt


ldm = c(paste0(lineage_defining_mutations$nt_alt,
                                      lineage_defining_mutations$nt_pos,
                                      lineage_defining_mutations$nt_ref),
                               paste0(lineage_defining_mutations$nt_ref,
                                      lineage_defining_mutations$nt_pos,
                                      lineage_defining_mutations$nt_alt))
        
minor_consensus_plotdata = patient_var_30 %>% #left_join(consensus) %>% 
    mutate(ref_nt = str_sub(ref_codon, 
                            start = codon_pos, end = codon_pos),
      ref_mutation = paste0(ref_nt, ntpos, minor),
      mutation = paste0(major, ntpos, minor), consensus = 
             ifelse(ntpos %in% lineage_defining_mutations$nt_pos, 
                    ifelse((mutation %in% ldm) | (ref_mutation %in% ldm),
                           TRUE, FALSE),      
                           FALSE)) %>% #filter(consensus == FALSE) %>%
  group_by(ntpos, consensus) %>% 
  summarize(consensus = consensus, n = n()) %>%
  unique %>% drop_na %>% ungroup %>% arrange(-n) %>%
  mutate(rank = 1:nrow(.))

`summarise()` has grouped output by 'ntpos', 'consensus'. You can override using the `.groups` argument.

minor_prevalence_across_genome_unlog <- minor_consensus_plotdata  %>%
    ggplot(aes(x = ntpos, y = n, fill = consensus, color = consensus)) + 
    # label your data too for the SALMON
    geom_bar(stat = "identity") + 
  scale_color_manual(values = c("salmon","grey")) + scale_fill_manual(values = c("salmon","grey")) +
    
      geom_point( aes(x = ntpos, y= n), shape = "") +
    #scale_color_manual(values = c("salmon","black")) + 
  theme_bw() + 
  geom_hline(yintercept = nrow(patient_counts_30)*0.01, linetype="dotted") + 
  xlab("Nucleotide position") + ylab("No. samples w/minor variant at site") + 
  annotate("rect", xmin=27894, xmax=28295, ymin=0, ymax=Inf, alpha=0.2, fill="#85D4E3") + 
  annotate("rect", xmin=28274, xmax=29533, ymin=0, ymax=Inf, alpha=0.2, fill="#F4B5BD") + 
  annotate("rect", xmin=21563, xmax=25384, ymin=0, ymax=Inf, alpha=0.2, fill="#FAD77B") +  xlim(0,29903)

minor_prevalence_across_genome_unlog

minor_prevalence_across_genome_unlog_labeled = 
  minor_prevalence_across_genome_unlog + 
  geom_text_repel(max.overlaps = Inf, ylim  = c(500,NA), size = 3, angle = 90, 
                  segment.linetype = 3, segment.size = 0.6, 
                  segment.curvature = -1e-20, 
                  arrow = arrow(length = unit(0.015, "npc"), 
                                type = "closed"),
                  data = minor_consensus_plotdata %>%
                    filter(rank < 50),
                    aes(x = ntpos, y = n, label = ntpos, color = consensus, fill = consensus))

Warning: Ignoring unknown aesthetics: fill

minor_prevalence_across_genome_unlog_labeled

genes <- fread("ntpos_gene_update.csv", data.table = F)

gene_limits<- genes %>% group_by(gene_id) %>% summarise(start=min(ntpos), end=max(ntpos)) %>% filter(gene_id!="") %>% mutate(molecule="")

gene_arrows<-ggplot(gene_limits, aes(xmin = start, xmax = end, y=molecule, label = gene_id)) +
  geom_gene_arrow(arrowhead_height = unit(3, "mm"), arrowhead_width = unit(1, "mm")) + geom_gene_label() + geom_segment(aes(x=266,y=1.5,xend=13468,yend=1.5), size=0.2) + xlim(0,29903) +
  annotate("text", label="ORF1a", x=5000, y=1.41, size=3) + geom_segment(aes(x=13468,y=1.4,xend=21555,yend=1.4), size=0.2) + annotate("text", label="ORF1b", x=18000, y=1.31, size=3) +
  theme_genes() + ylab(NULL) + xlab(NULL) 

Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
Please use `linewidth` instead.

((
  genome_fig_unlog<-plot_grid(gene_arrows, minor_prevalence_across_genome_unlog_labeled + theme(legend.position = "bottom"), nrow=2, rel_heights=c(0.25,1), axis="lr", align="hv")
))



ggsave("ggsave/genome_fig_unlog.pdf", plot = genome_fig_unlog, height = 5, width = 12)

Characteristics of highly recurrent minor variants

# What's the range of minor allele frequencies the highly recurrent minor variants are found at?
highly_shared_sites_top_ranked = highly_shared_sites_ranked %>% filter(rank <= 35)
highly_shared_sites = patient_var_30_for_rank %>% filter(paste0(gene,".",label) %in%
                                     paste0(highly_shared_sites_top_ranked$gene, ".", 
                                            highly_shared_sites_top_ranked$label)) %>%
          group_by(label, ntpos) %>% 
  summarize(ntpos, ntpos_count = n()) %>% distinct() %>% 
  arrange(-ntpos_count) %>% 
  filter(ntpos_count > (patient_var_30$MCoVNumber %>% unique %>% length)*0.01) %>% pull(ntpos)

`summarise()` has grouped output by 'label', 'ntpos'. You can override using the `.groups` argument.

((
  hss_maf<-patient_var_30 %>% filter(ntpos %in% highly_shared_sites) %>% 
    ggplot(aes(x = as.factor(ntpos), y = minorfreq)) + 
    geom_point(alpha = 0.2, position = position_jitter(width = 0.1)) + theme_bw() + 
    geom_boxplot(outlier.shape = NA, alpha =0.5) + 
    theme(axis.text.x = element_text(angle=90)) + #geom_boxplot(outlier.alpha = 0) + 
    xlab("Nucleotide position") + ylab("Minor allele freq") + 
    scale_y_continuous(trans="log2", breaks = c(0.01*2^(0:5), 0.5))
  ))

ggsave("ggsave/hss_maf.pdf", width = 4, height = 1, plot = hss_maf)

((
  hss_maf_unlog<-patient_var_30 %>% filter(ntpos %in% highly_shared_sites) %>% 
    ggplot(aes(x = as.factor(ntpos), y = minorfreq)) + 
    geom_point(alpha = 0.2, position = position_jitter(width = 0.1)) + theme_bw() + 
    geom_boxplot(outlier.shape = NA, alpha = 0.5) + 
    theme(axis.text.x = element_text(angle=90)) + #geom_boxplot(outlier.alpha = 0) + 
    xlab("Nucleotide position") + ylab("Minor allele freq") #+ 
    #scale_y_continuous(trans="log2", breaks = c(0.01*2^(0:5), 0.5))
  ))
ggsave("ggsave/hss_maf_unlog.pdf", width = 4, height = 2, plot = hss_maf_unlog)

#How many different changes do you find at each recurrently-mutated position? How many such mutations are a reversion to the reference?
hss_info <- fread('hss_data2.csv', data.table = F)

((
  hss_plot = patient_var_30_for_rank %>% mutate(label_gene = paste0(gene, ": ", label)) %>%
    filter(ntpos %in% highly_shared_sites) %>% 
    group_by(ntpos, major, minor) %>% mutate(n=n(), median_frequency = median(minorfreq)) %>%     mutate(change=paste0(major,">",minor)) %>% 
    mutate(ref_reversion=if_else(minor==refnt,"yes","no")) %>% 
    mutate(median_MAF=if_else(median_frequency<0.02, 
                              "low (median <2% MAF)","high (median >=2% MAF)")) %>%
    select(change, n, ref_reversion, ntpos, label_gene, median_MAF) %>% distinct %>%
    ggplot(aes(x=change, y=n, color=ref_reversion)) + 
    geom_bar(stat="identity", aes(fill=median_MAF)) + 
    scale_color_manual(values=c("white","red")) + 
    scale_fill_manual(values=c("darkmagenta","steelblue")) + 
    facet_wrap(ntpos~label_gene, scales="free") + 
    theme_pubr() + theme(axis.text.x=element_text(angle=90))
))

Adding missing grouping variables: `major`, `minor`

ggsave('ggsave/hss_plot.pdf', plot = hss_plot, width = 11, height = 13)

ldm_vector = lineage_defining_mutations$nt_pos
heatmap_spectra_tmp = patient_var_30_for_rank %>% 
  mutate(label_gene = paste0(gene, ": ", label)) %>%
    filter(ntpos %in% highly_shared_sites) %>% 
    group_by(ntpos, major, minor) %>% mutate(n=n(), median_frequency = median(minorfreq)) %>%
  mutate(change=paste0(major,">",minor), 
         transition = ifelse(change %in% c("C>T","T>C","A>G","G>A"), 1, 0)) %>%
     mutate(ref_reversion=if_else(minor==refnt,1,0)) %>% ungroup() %>%
  select(change, n, ntpos, label_gene) %>% distinct %>% 
  group_by(ntpos) %>% mutate(ntpos_n = sum(n)) %>% ungroup() %>%
  mutate(n_prop = n/ntpos_n) %>% select(-n) %>% 
  mutate(ldm = as.factor(ifelse(ntpos %in% ldm_vector, 1, 0))) %>%
  mutate(label = paste0(label_gene, " - ", ntpos, " (", ntpos_n, ")"))


annotation_row_df = heatmap_spectra_tmp %>% select(label, ldm) %>% 
  distinct %>% as.data.frame() %>%
  column_to_rownames("label")
annotation_row_df

#%>%
heatmap_spectra = heatmap_spectra_tmp %>% select(-ldm) %>% spread(change, n_prop) %>% 
  arrange(-ntpos_n) %>% distinct() %>%  
  select(-c(ntpos, label_gene, ntpos_n)) %>% replace(.,is.na(.),0) %>%
  column_to_rownames("label")

annotation_col_df = heatmap_spectra_tmp %>% select(change) %>% distinct %>% 
  mutate(transversion = 
           as.factor(ifelse(!(change %in% c("C>T","T>C","A>G","G>A")), 
                            1, 0))) %>% 
  column_to_rownames("change")

maf = heatmap_spectra_tmp %>% select(ntpos, label) %>% distinct %>%
  left_join(patient_var_30 %>% select(ntpos, minorfreq)) %>% 
  select(label, minorfreq)

Joining, by = "ntpos"

maf_list = split(maf$minorfreq,maf$label)

row_ha = rowAnnotation(df = data.frame(
  ldm = annotation_row_df[rownames(heatmap_spectra),]), 
  col = list(ldm = c(`1`="grey", `0`="white")),
  MAF = anno_boxplot(maf_list, pch = 20, size = unit(0.5, "mm")))


heatmap_spectra_reorder = heatmap_spectra %>% select(c("A>G","G>A","C>T","T>C",
                                               "A>C", "C>A", "A>T", "T>A",
                                               "C>G", "G>C", "G>T", "T>G"))
col_ha = columnAnnotation(df = data.frame(
  transversion = annotation_col_df[colnames(heatmap_spectra_reorder),]), 
                          col = list(transversion = c(`1`="grey", `0`="white")))

heatmap_spectra_hss = Heatmap(heatmap_spectra_reorder,
                              name = "prop/site",
                              col = colorRampPalette(
                                c("white", "brown4"))(100),
                              show_row_dend = F, 
                              cluster_columns = F, left_annotation = row_ha,
                              rect_gp = gpar(col = "grey90", lwd = 1),
                              top_annotation = col_ha,
                              heatmap_legend_param = 
                                list(direction = "vertical")) 

Warning: The input is a data frame-like object, convert it to a matrix.Warning: Values in row annotation 'MAF' have a different order of names from the matrix
row names. It may lead to wrong conclusion of your data. Please double check.

pdf(qq("ggsave/heatmap_spectra_hss.pdf"), width = 6.5, height = 7)
draw(heatmap_spectra_hss, merge_legend = TRUE, 
     heatmap_legend_side = "right", 
    annotation_legend_side = "right")
dev.off()

null device 
          1 

heatmap_spectra_hss

# 
#           border_color = "grey90", colorRampPalette(c("white", "brown4"))(100),
#          treeheight_row = 0, treeheight_col = 10, annotation_col = annotation_col_df,
#          annotation_row = annotation_row_df, annotation_colors = annotation_colors)
# pacman::p_load(ComplexHeatmap)
         

major_minor_table(consensus_runs)

`summarise()` has grouped output by 'run'. You can override using the `.groups` argument.`summarise()` has grouped output by 'run'. You can override using the `.groups` argument.Joining, by = c("run", "ntpos")

null device 
          1 

major_minor_table(consensus_runs)

null device 
          1 

major_minor_table(consensus_runs, all = T, 
                  outfile = "ggsave/heatmap_run_mutation_all.pdf",
                  threshold_sites = minor_consensus_plotdata$ntpos[1:1000] %>% unique)

null device 
          1 

#all minor variants that are found in at least 2% of high-coverage samples in Lythgoe et al https://www.science.org/doi/suppl/10.1126/science.abg0821/suppl_file/abg0821-lythgoe-sm.pdf
lythgoe_hss<-c(29320, 25628, 25807, 20364, 357, 25627, 25532, 11083, 239, 238, 356, 16740, 20989, 19393, 15009, 21826, 19937, 26289,28253,25507,29862,11052,29864,25529,369,22453,13571,29538,21575,75,20993,19473,15474,635,1226,29747,24933,19210,13303,28936,26334,16887)

tonkin_hss<-c(11075L, 11083L, 11074L, 522L, 26780L, 1547L, 14408L, 28253L, 13914L, 26137L, 635L, 241L, 26785L, 683L, 23403L, 28881L, 13778L, 25521L, 14805L, 9491L, 11071L, 13780L, 21575L, 29242L, 9430L, 203L, 3037L, 686L, 3096L, 16375L, 16466L, 17167L, 23559L, 9438L, 21137L, 29051L, 23555L, 26333L, 9474L, 13571L, 16887L, 26787L, 28603L, 24213L, 11651L, 27925L, 514L, 26781L, 558L, 1912L, 9479L, 11073L, 12578L, 26144L, 26681L, 29241L)
intersect(lythgoe_hss, highly_shared_sites)
intersect(tonkin_hss, highly_shared_sites)

homoplasic_site_list<-c(187L, 1059L, 2094L, 3037L, 3130L, 6990L, 8022L, 10323L, 10741L, 11074L, 11083L, 13408L, 14786L, 15324L, 19684L, 20148L, 21137L, 21575L, 24034L, 24378L, 25563L, 26144L, 26461L, 26681L, 28077L, 28826L, 28854L, 29700L)
intersect(homoplasic_site_list, highly_shared_sites)

tally_in_minors<-patient_var_30 %>% filter(ntpos %in% highly_shared_sites) %>% group_by(ntpos) %>% summarise(n_minor_variants_houston=n()) %>% mutate(freq_minor_variants_houston=n_minor_variants_houston/nrow(patient_counts_30))

tally_in_gisaid<-fread('hss_gisaid_tallies.csv', data.table = F) %>% select(ntpos, n_consensus_variants_gisaid=n) %>% mutate(freq_consensus_variants_gisaid=n_consensus_variants_gisaid/3109421) #USA gisaid sequences as of 6/13/2022
#SUPP FIG 10

tally_in_minors_gisaid = fread("processing/tally_in_minors.txt", data.table = F) %>% 
  left_join(minor_consensus_plotdata) %>% mutate(ldm = as.factor(consensus))

tally_in_minors_gisaid
((
  gisaid_vs_hmh_plot = #tally_in_minors %>% left_join( tally_in_gisaid) %>% 
    
    tally_in_minors_gisaid %>%
    ggplot(aes(x=freq_minor_variants_houston, y=freq_consensus_variants_gisaid, 
               color = ldm)) + geom_point(alpha=0.5) + 
    theme_pubr(legend= "right") + 
  xlab("Fraction of samples with minor variant in Houston data") + 
    ylab("Fraction of samples with \n consensus variant in GISAID Global") + 
    geom_text_repel(aes(label=ntpos), max.overlaps = 10, label.size=0.1) + 
    geom_abline(slope = 1, intercept = 0, linetype="dashed", color = "gray") +
    scale_y_continuous(trans="log10") + scale_x_continuous(trans="log10") + theme(legend.position = "bottom")
))

gisaid_vs_hmh_marginal = ggMarginal(gisaid_vs_hmh_plot,
             groupColour = T,
  groupFill = T)

ggsave("ggsave/gisaid_vs_hmh_plot_marginal.pdf", plot = gisaid_vs_hmh_marginal, width = 6, height = 5)

# MUTATION SPECTRA ANALYSIS

# reproducibility analysis from RMD 01
plot_nonreproducible_spectra_heatmap = readRDS(file = "plot_nonreproducible_spectra_heatmap.rds")
plot_reproducible_spectra_heatmap = readRDS(file = "plot_reproducible_spectra_heatmap.rds")

hmap<-table(patient_var_30$major[patient_var_30$major!=""], patient_var_30$minor[patient_var_30$minor!=""]) %>% data.frame() %>%
  ggplot(aes(x=Var1, y=Var2, fill=Freq/sum(Freq))) + geom_tile(colour = "black") +
  scale_fill_gradient(low = "white",
                      high = "darkblue") +
  theme_minimal() +   labs(fill = "Number",
       x = "Consensus allele",
       y = "Minor allele", title = "Post-sample filter (Ct, n_var)") #Most frequently C>T mutations

# Now beecause some peaks are really prevalent like the twin peaks at 29187/8, 
# it's important to not factor that in too much. So let's just assign those mutations
# to equate to the median number a mutation appears which is 1.

mutation_spectra = patient_var_30 %>% unite(col = "spectra_mutation", 
                                            major, ntpos, minor, remove = FALSE)
mutation_spectra_counts = mutation_spectra %>% count(spectra_mutation)

mutation_spectra_counts %>% ggplot(aes(y=n)) + geom_boxplot() + 
  scale_y_continuous(trans = "log1p")

mutation_spectra_counts$n %>% quantile() # median is 1

mutation_spectra_unique = mutation_spectra %>% 
  select(major, minor, spectra_mutation) %>% unique()

hmap_unique <-table(mutation_spectra_unique$major[mutation_spectra_unique$major!=""], mutation_spectra_unique$minor[mutation_spectra_unique$minor!=""]) %>% data.frame() %>%
  ggplot(aes(x=Var1, y=Var2, fill=Freq/sum(Freq))) + geom_tile(colour = "black") + # grid color
  scale_fill_gradient(limits = c(0,0.4), low = "white",
                      high = "darkblue") +
  theme_minimal() +   labs(fill = "Fraction",
       x = "Consensus allele",
       y = "Minor allele", title="Unique minor variants")

rescale_fill = function(hmap_unique) {
  hmap_unique_rescaled = hmap_unique + 
    scale_fill_gradient(limits = c(0,0.4), low = "white",
                      high = "darkblue") +
    geom_label(fill = "white", alpha = 0.3, aes(x = Var1, y = Var2, label = round(Freq/sum(Freq), digits = 2)))
  return(hmap_unique_rescaled)
}

((
  rep_nucleotides = ggarrange(rescale_fill(plot_nonreproducible_spectra_heatmap), 
                              rescale_fill(plot_reproducible_spectra_heatmap), 
                              rescale_fill(hmap), 
                              rescale_fill(hmap_unique), common.legend = T)
))
rep_nucleotides

ggsave("ggsave/heatmap_spectra_replicate_variants.pdf", rep_nucleotides, height = 6, width = 6)